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Abstract 

In earlier work we have studied a method for discretization in 
time of a parabohc problem which consists in representing the ex- 
act solution as an integral in the complex plane and then applying 
a quadrature formula to this integral. In application to a spatially 
semidiscrete finite element version of the parabolic problem, at each 
quadrature point one then needs to solve a linear algebraic system 
having a positive definite matrix with a complex shift, and in this pa- 
per we study iterative methods for such systems. We first consider the 
basic and a preconditioned version of the Richardson algorithm, and 
then a conjugate gradient method as well as a preconditioned version 
thereof. 
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1 Introduction 

Let be a complex finite-dimensional inner product space, and let A be a 
positive definite Hermitian linear operator in V, with spectrum We 
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shall consider iterative methods for the linear equation 



zw + Aw = g, where z = x + iy ^ —a{A). 



(1.1) 



Such equations, with a complex shift z of the positive definite operator A, 
need to be solved in a method for discretization in time of parabolic equations, 
based on Laplace transformation and quadrature, which has been studied 
recently, as will be made more specific below. Equations of the form (1.1) 
arise also from the spatial discretization of the Helmholtz equation, cf. [6], 
but in that context the ^-vahics typically of interest differ from those we wish 
to consider; for our application to the heat equation, arg^ is bounded away 
from ±7r. In this paper we shall consider a basic Richardson iteration and a 
conjugate gradient (CG) method for (1.1), as well as preconditioned versions 
of these methods. Another approach, not discussed here, is to reformulate the 
complex linear system as an equivalent real one with twice as many equations 
and unknowns, cf., e.g., [2] and the list of references therein. 

We begin by sketching the time discretization method referred to above. 
In a complex Hilbert space H we consider the initial-value problem 



where ^4 is a positive definite Hermitian operator in 7^. To represent its 
solution we apply the Laplace transform, writing 



Ut + Au — f{t), for i > 0, with ii(0) = Wo, 



(1.2) 




Under appropriate assumptions on f{t) we then have formally 



{zl + A)w{z) = Wo + f{z) 



or, multiplying by the resolvent of —A, 



w{z) = R{z)g{z), where R{z) 



{zl + A) 
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Applying the inverse Laplace transform, we obtain 




for := {z : Rez = u} and u > 0. With ip e (^tt, tt) and T a new contour 
in = : | a.Ygz\ < (p}, homotopic with F^,, we may write 
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A suitable parametrization of F, written z — z(^) for ^ G R, yields 

/oo 1 
v{i,t)di, where v{i,t) _ e^«)*^(^(0)/(0- (1-3) 

We assume Rc 2;(^) — >■ — oo as |^| — )■ oo so that e^^^^* — )■ 0, for i > 0. 

We now define an approximate solution of (1.2) by means of an equal- 
weight quadrature rule, applied to the integral in (1.3), 

V,{t) := h v{^„ ^) = ^ E ^'Mzj)z'j, (1.4) 
j=-q j=-<i 

where, for an appropriate k > 0, we have set 

^j-.^jkeR, Zj:^z{Q, z'j-.^z'iQ, for|i|<g. (1.5) 

To compute Ug{t), we need to solve the 2q + 1 "elliptic" equations 

(zjl + A)w(zj) = g(zj), for \j\ < q. 

These equations are independent, and may thus be solved in parallel. We note 
that the w{zj) determine Uq{t) for all t > 0, but we can expect an accurate 
approximation only for t in some restricted interval that depends on the 
choice of the quadrature step k and of the parametric representation z{^). 

In our presentation we shall follow the analysis of [11]. Specifically, we 
use for r the left branch of the hyperbola (a; — 1)^ — = 1 in the complex 
plane, parametrized by 

z{^) = 1 - cosh ^ + i sinh ^, ^ G M, 

and take k — log q/q for the step size in (1.4). This means that 



Zj — Xj + iyj — 1 — cosh ^iZi^^^^ -\- i sinh ^iZi^^I^^ ^ for \j\ < q. 

In particular, Zg — 1 — {q + q~^)/2 + i {q ~ 9^^)/2 ~ — q'/2 + iq/2 for large q. 
Under the appropriate assumptions about the data of the problem we then 
have the error estimate, see [11], with < tg < T < oo. 

We now want to apply this time discretization scheme to the semidiscrete 
finite element approximation of the heat equation, with elliptic operator Lu — 
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— V • (aVit), and consider thus the initial boundary- value problem for u — 

u(x,t), 

Ut + Lu = /(-, t), in fl, with u = on dfl, for t > 0, 
u{(J, ■) = Mo, m \l, 

where, for simplicity, we will assume that the diffusivity a is a (positive) 
constant, and that f2 is a convex polygonal domain in M^. This problem 
is the special case of (1.2) with H = L2{fl) and A = L, taking D{L) = 

Let {Vh} C Hq{Q) be a family of piecewise linear finite element spaces, 
based on a family of regular triangulations Th = {r} of fl. With (v, w) — 
J^vwdx, the standard Galerkin, spatially semidiscrete approximation of 
(1.6) is 

{uh,u X) + a^^Uh: Vx) = (/, x), Vx e T4, t > 0, with Uh{0) = Uoh, 

where, with Ph : L2{fl) Vh the L2-projection onto Vh, we may take, e.g., 
uoh — PhUQ- Introducing the discrete elliptic operator Lh Vh ^ Vh, defined 

by 

(L^V',X) =«(VV',Vx), yi^,x&Vh, 

the spatially semidiscrete initial-value problem may also be written 
Uh,t + LhUh = Phfi-, t), for t > 0, with ^,^(0) = P^Mq, 

which is of the form (1.2) with = Vh, equipped with the L2 inner prod- 
uct, and A = Lh- The fully discrete solution defined by our above time 
discretization method (1.4) now takes the form 




j=-Q 



with Zj, Zj as in (1.5) and where the Wh{zj) are derived from 
{zjl + Lh)wh{zj) = Phg{zj), for \j\ < q, 

or, in weak form, 

Zj{wh(zj),x) +a{Vwh(zj),Vx) = {g(zj),x), V x e K- (1-8) 

As before, these problems may be solved in parallel. We note that they 
are special cases of (1.1), with V = Vh and A — Lh- Under appropriate 
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assumptions on the data [11] the error in the fully discrete solution may be 
bounded as 

\\U,,,{t)-u{t)\\ < Ct^A^o^rnh' + e-'^/'^^"), for t G [to,T] C (0,oo). 

(1.9) 

To express (1.8) in matrix form, let {Pi}^i be the interior nodes of Th 
and the associated nodal basis functions, so that v & Vh may be 

written as v — ^iLi'^i^i with Vi v{Pi). Let Ai = {ran) and S — 
(sii) be the mass and stiffness matrices, where mu := ($j, $;) and sn :— 
a(V$i, V$;), respectively. With w — Wh{zj) — X^^i Wi^i, equation (1.8) is 
then equivalent to 

ZjMw + Sw — g or Zj w + M~^Sw = M~^g, (1.10) 

where the components of the load vector arc = ((yf(zj), $«) . The second 
equation in (1.10) is of the form (1.1) with Av = Ai~^Sv and g = Ai~^g. 
However, instead of the standard unitary inner product {v, w) = J2iLi '^i'^u 
we equip V = with (v, w) — {M.V, w) so that A is Hermitian: {Av, w) — 
{AiAv,w) = {Sv,w). In our study of iterative methods for (1.10), we 
develop the theory for an abstract operator A satisfying our assumptions, 
and discuss separately the practical implications for the specific choices A — 
Lfi and, especially, A — M.~^S. 

As an alternative to the standard Galerkin method we may consider the 
lumped mass modification, in which the mass matrix A4 is replaced by a 
diagonal matrix V; we refer to [18] for details. 

For any Hermitian operators A and B in V, with B positive definite, we 
will write Xj = Xj{A, B) for the jth generalized eigenvalue of A with respect 
to B, that is, Avj — XjBvj with Vj ^ 0. We order these eigenvalues so that 
Ai < A2 < ■ ■ ■ < Aat, and use the abbreviation Xj{A) = Xj{A, I) 

The program for time discretization of parabolic equations sketched above 
was initiated in Sheen, Sloan and Thomee [15, 16], and continued in Gavri- 
lyuk and Makarov [7], McLean, Sloan, and Thomee [11] and McLean and 
Thomee [12, 14, 13], cf. also Thomee [17], and the error in (1.7) was ana- 
lyzed in both 1/2(0) and Loo(f^), under various assumptions on the data of 
the problem. In the latter papers also fractional order diffusion equations 
were treated. 

In these papers the analysis was illustrated by numerical examples. These 
were carried out in simple cases, in one space dimension and also in the case 
of a square spatial domain in two dimensions, and direct solvers were used 
for the linear system (1.8). However, even though powerful direct solvers are 
available, for large size problems in more complicated geometries, particularly 
in 3D, it may be natural to apply iterative methods, and our purpose in this 
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paper is therefore to begin a study of such methods for equations of the 
form (1.1), with apphcation to the heat equation in mind. Some preliminary 
results on this problem were sketched in [16], using the Richardson iteration 
algorithm for (1.8) and for a preconditioned form of this equation, and in 
Section 2 below we extend and improve these results. 

Prom a knowledge of the extremal eigenvalues of ^4, we can determine the 
optimal value of the complex acceleration parameter, optimal in the sense of 
minimizing the error reduction factor of the Richardson iteration. For the 
finite element problem on quasiuniform triangulations with maximal mesh- 
size h, the basic Richardson method converges slowly, with the error in the 
nth iterate bounded by (1 — c/i^)", for c > depending on z, but with the 
convergence rate improving with growing \z\. We also study preconditioned 
versions of this method, first using the special preconditioner Bz = {^zl + 
A)~^, where ^z > — Ai(A), which may be analyzed in the same way as 
the basic method, and we show that the error reduction factor is bounded 
away from 1 as Ajv — > oo. We then consider a general preconditioner Bz, 
and prove geometric convergence in the norm \[v]\ :— {B~^v,vy^^, where 
the acceleration parameter is defined in terms of bounds for the spectrum 
of Bzii^zl + A). 

In Section 3 we analyze a CG method, which does not involve choosing 
an acceleration parameter. Generalizing the usual convergence analysis to 
allow the complex shift of A in (1.1), we show geometric convergence of the 
iterates Wn, which follows from the error bound 

sec(iarg2;) Xi + Xn + 2z 
\\\w„ - w\\\ < ,„ . |||wo - w|||, with := — , (1.11) 

\-ln[Sz)\ Aat - Ai 

where |||f||P := + {Av,v) and T„ is the Tchebyshev polynomial of 

degree n, and where Xj = Xj{A). Since T„(s^) = ^{rj^ + ?7~") with \r]z\ < 1, 
this indicates geometric convergence with rate \r)z\"'. For the finite element 
problem discussed above, we find that jr^^l < 1 — ch with c > depending 
on z, giving a better convergence rate than Richardson iteration. 

If the equation is preconditioned with Bz = {^zl + A)~^ , for appropri- 
ate iiz-i and if we let z = {z — Hz)~^, then the preconditioned equation is 
equivalent to zw + BzW — zBzQ, which again has the form (1.1), and a sim- 
ilar convergence result holds with an error reduction factor bounded away 
from 1 as Aat — )■ oo. 

It is natural to consider more general preconditioncrs also for the CG 
iteration. The preconditioned equation zBzW + BzAw = BzQ is again equiv- 
alent to an equation of the form (1.1), namely, zv -\- bI^^ABz^^^v — Bl^^g, 
where v — bI^^w and the transformed operator bI is Hermitian 



6 



and positive-definite with respect to the inner produce [v,w] — {B^^v,w). 

4-1/2 

However, computing the action of Bz will usually be costly, so we instead 
work with the preconditioned equation in its original form. Although the 
error is still optimal in a certain sense, we are not able to show a precise 
error bound of the type (l-H). 

Section 4 develops the algorithmic implementation of the CG method. 
For the basic method, the successive iterates satisfy a three term recursion 
relation. The same is true of the preconditioned method for the special 
choice Bz — (//^/ + A)~^, but not necessarily for a more general precondi- 
tioner. 

Use of an iterative solver means that we compute an approximation Wh{zj) 
in place of the true finite element solution Wh{zj), so that in place of (1.7) 
we obtain 




j=-9 



If \\wh{zj) - Wh{zj)\\ < Ej, then 

E(t) \\U,,n{t) - U,,,(t)\\ < 1^ E 141' (1-12) 

j=-<i 

and we may use this estimate as the basis for a stopping criterion. In view 

of the error estimate (1.9) we see that it is desirable to choose the solver 
tolerance Ej in such a way that S{t) < C{h^ + 6^^^^^°^''). The presence of the 
factor e^^^\z'j\ allows £j to increase with |j|; see (5.1) below and remember 
that < 0. 

In the final Section 5 we illustrate our error analysis by numerical calcu- 
lations in a concrete case of (1.6), and discuss how to choose the parameters 
to balance the contributions to the error of the discretizations in space and 
time and in the iterative procedure. 

2 Iteration algorithms of Richardson type 

We now assume, as in (1.1), that ^4 is a positive definite Hermitian operator 
in a finite-dimensional complex inner product space V . with extremal eigen- 
values Ai = \i{A) and = Ajv(A), and for brevity put Az := zl + A. In 
this section, following [16], we consider first the basic Richardson iteration 
with acceleration parameter a e C, applied to AzW = g, 

^ [I - aAz)w^ + ag, for n > 0, with w° given. (2.1) 
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The error reduction in each time step is then described by the inequahty 

- w\\ < \\I - aA^W - wll, 
and since is a normal operator in V, 

\\I - aAJ\ ^ max \1 - a(z + \)\. (2.2) 

In (2.1), in addition to choosing the issue is to select a e C so that 
the norm in (2.2) is as small as possible. For ^ = 0, as is well known, the 
optimal choice of a is 2/(Ai + Xn), which gives 

k(A) - 1 , , ^, Aat 
K[A) + 1 Ai 

When A — Lfi is based on a quasi-uniform family of triangulations 7h we 
have K,{A) — 0{Xn) — 0{h~'^) and hence, in this case, 

||/ - a^ll < 1 - c/^^ with c> 0. (2.3) 

To determine an optimal a in (2.2), we shall have use for the following 
lemma. 

Lemma 2.1. Zet a, b G C 6e nonproportional, and [a, b] C C the line segment 
with endpoints a and b. Set 

F(a) :— max 11 — aX\, where a & C 

AG[a,b] ' 

Then F{a) < 1 for suitable a, and F{a) is minimized by 

a — — ^ — r, where c := ^(a + b), d := i(h — a), 
c + sd 2v " V j> 

and where s e M minimizes the real rational function 

a 2 jdps^ ^2sRe((c- a)d) + |c- ap 



Ris) 



1 - 



c + sd 



|d|V + 2sRe(cd) + |c|' 



The minimizing value of s is given, with the ± sign being that o/Re(ad), by 



where _ 

2Re(ac)-|aP , , 2fi Re(cd) - IcI 

fi := = and U '■= tttt; 

2Re(ad) |d|2 
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Proof. We first note that a may be chosen so that F{a) < 1. In fact, we 
may first rotate the hne segment [a, b] around the origin so that it becomes 
parallel to and to the right of the imaginary axis, which determines argci, 
and then shrink the line segment thus rotated so that it comes inside the 
disk 1^; — 1| < 1, giving |a|. 

For a to be optimal, we must have |1 — aaj = |1 — abj, and thus also 
|1/q; — a| = |1/q; — b|. Therefore, 1/a has to be chosen on the line in C 
through the midpoint c = |(a + b), which is perpendicular to b — a, or has 
the direction of d = i{b — a), so that 1/a = c + sd, or a = l/(c + sd), 
with s G M, and R{s) = |F(a(s)) 

Since d 7^ we have R{s) ^ 1 as s — )■ ±00, and if Re(ad) > ( < 0, 
respectively) then R{s) < 1 ( > 1, respectively) for large s > 0. Note that 
since a = ai + ia2 and b = bi + ih2 are nonproportional, we have Re(ad) = 
— Re(ia(b — a)) = — Re(iab) = a2bi — aib2 7^ 0. A simple calculation shows 



R'{s) = 



2Re(ad)|d| 




(S2| 


|d| 


|2 + 2sRe(cd) + 


c 


2^2 



SO that R{s) has just one maximum and one minimum, with the maximum 
to the left of the minimum if and only if Re (ad) > 0. □ 

We are now ready to show the following. 

Theorem 2.2. Let z — x + iy with axgz e (— 7r,7r) and determine a — by 
taking a — z + Xi and b = z + in Lemma 2.1. Then, for Ajv sufficiently 
large, the error reduction factor in (2.1) satisfies 

Ez :— min ||7 — aA^W = |1 — a^a] < 1 — cX]^^, with c — c{z, Ai) > 0. 

a 

Proof. With the notation of Lemma 2.1 we have, for A^v — > 00, 

2Re(ac) = (a;+Ai)Ajv+0(l), Re(ad) = yXN+0{l), Re(cd) = yA^+0(l). 

Hence, putting s± := - ( - a; - Ai ± ■\/{x + Ai)^ + t/2) / (^2y), 

h = ^^ + 0{X],'), f, = -^ + 0{X],'), Smin = s± + 0(A]^^), 

and it follows by Lemma 2.1 that ct^ = (| + is±)~^X]^^ + 0{X]^^). Therefore, 
|1 - a,a\'^ = 1-2 Re(Q;^a) + \a,a\'^ = 1 - ^A]^^ + 0{X]^^), 
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where 



V 2 + ^'^± / 



^ x + Ai + 2ys± ^ ^J{x + \if + y^ 

since the sign in ± is that of y, and the desired estimate follows for Aat 
sufficiently large. □ 

When, as above A = Lh, with {Th} quasiuniform, so that \n ^ ch^^, the 
error bound is of the same form as in (2.3), except that now the constant c 
depends on z. 

The rate of convergence shown in Theorem 2.2 is too slow for the iteration 
to be of practical use. In Table 2.1, we show the values of the parameter a — 
pe~^'^ and the error reduction factor given by Theorem 2.2, with z = Zj 
on the hyperbola {x — 1)^ — = 1, for even j in the range < j < g = 20. 
Here, the operator A is from the model problem described in Section 5, for 
which Ai 1 and Aat ~ 4, 000. 

One way to improve the convergence of the iterative method (2.1), con- 
sidered briefly in [16], is to precondition the linear system by multiplication 
by a positive definite Hcrmitian operator B^, which, in contrast to the choice 
in [16], we here allow to depend on z. Rewriting (1.1) as 

GzW = := B^g, where := B^A^, (2.4) 

the Richardson iteration algorithm becomes 

^n+i ^ (J _ aG,)w'' + ag,. (2.5) 

We first consider the special preconditioncr B^ = {fizl + where 
/iz > — Ai. One could choose, for example, /x^ = 0, as in [16], or //^ = \z\. 
For /iz — we have B^ — A~^, independently of z, and ior /i^ — \z\, Gz is 
bounded in z. Since 

Gz = Gz{A, i^z) = {l^zl + A)-\zI + A), (2.6) 

the error reduction is now measured by 

z -\- \ 

\\I - aGz{A,ij,z)\\ = max \1 - aGz{X,iiz)\, Gz{X,Hz) = — T' (2-7) 

AecT(A) + A 

and we want to choose a so that this quantity is as small as possible. 
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Theorem 2.3. Let z — x + iy with aigz e (— 7r,7r), let /iz > — Ai, and de- 
termine a — az by taking a = Gz{\i,iiz) and b = Gz{\n-iIJ'z) Lemma 2.1. 
Then the error reduction factor in (2.5) is bounded independently of by 

Sz := ||/ - azGz{A,[iz)\^ = |1 - olzA < c{z,\x,[iz) < 1- (2.8) 
Proof. We note that 

Gz{X,iiz) e [Gz{Xi,ij,z),Gz{Xn,ij>z)], for A e [Ai, Ajv], 

and that Gz{Xn,IJ'z) ^ 1 as X^ oo. Thus, Gz{X,iJ,z) G [a, 1] for all 
A e o-{A), and since this is a fixed line segment, Lemma 2.1 shows the 
theorem. □ 

Since, for z, Ai and Xn given, the factor is an explicit, albeit compli- 
cated, function of /iz, it is natural to choose /iz as the value that minimizes 
this function. The numerical values of /x^ used in this section were determined 
in this way, via an optization routine, scipy . optimize . f minbound [9], based 
on a well-known algorithm due to Brent that docs not require derivative val- 
ues. We obtained almost identical results, not shown here, by setting 6=1, 
corresponding to Aat = oo. 

In Table 2.1, we see the dramatic effect of the preconditioner Bz — {iJ>zI + 
A)~^ on the error reduction factor, in the case of the model problem from 
Section 5, with z = Zj. Notice that Sz increases with j, whereas Sz decreases. 

Since computing the action of {nz^ + A)^^ is expensive, we now want 
to consider a more general preconditioner B^ (still assumed to be positive 
definite and Hermitian). Suppose first that z — and write B — Bq. If B~^ 
is spectrally equivalent to A, that is, if 

m{B-^v,v) <{Av,v) <M{B-^v,v), \/v eV, (2.9) 

for some positive m and M, then for suitable a the iterative scheme converges 
geometrically with respect to a suitable energy norm. More precisely, setting 

[vM-={B-\,w), \[v]\:=[vM"\ (2.10) 

the operator BA is Hermitian with respect to [■,■], with eigenvalues Xj{BA) = 
Xj{A,B~^) in the closed interval [m,M], so that k{BA) < M/m, and thus 

\[I -aBA]\ = <— if — ^^r—. 2.11 

^' k{BA) + 1- M + m Xx{BA) + Xn{BA) ^ ' 

In the general case of (2.4) with 2; 7^ 0, we shall write Gz in the form 

Gz^ BzAz = zBz + Bz{iizl + A), whexe z -.^ z - jiz- (2.12) 
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Table 2.1: Richardson iteration with a — pe and preconditioning with 
B, = {iiJ + A)-\ 





Theorem 2.2 


Theorem 2.3 


J 


T ■ 


7/ ■ 

yj 




Hz 


H^z 


c- 

'^Z 


n 

Hz 


H^z 


Hz 


c- 

'^Z 





0.00 


0.00 


4 


99e-04 


-0.00 


0.9995 


1.000 


-0.00 


0.00 


0.000 


2 


-0.05 


0.30 


4 


93e-04 


0.15 


0.9995 


0.988 


0.15 


0.00 


0.152 


4 


-0.18 


0.64 


4 


73e-04 


0.33 


0.9995 


0.947 


0.33 


0.03 


0.321 


6 


-0.43 


1.02 


4 


31e-04 


0.53 


0.9996 


0.864 


0.53 


0.16 


0.503 


8 


-0.81 


1.51 


3 


76e-04 


0.72 


0.9996 


0.753 


0.72 


0.51 


0.658 


10 


-1.35 


2.12 


3 


24e-04 


0.86 


0.9995 


0.650 


0.86 


1.14 


0.760 


12 


-2.10 


2.93 


2 


85^04 


0.96 


0.9995 


0.571 


0.96 


2.11 


0.821 


14 


-3.13 


4.01 


2 


58e-04 


1.03 


0.9994 


0.516 


1.03 


3.52 


0.856 


16 


-4.54 


5.45 


2 


39e-04 


1.07 


0.9993 


0.478 


1.07 


5.47 


0.878 


18 


-6.45 


7.38 


2 


25e-04 


1.10 


0.9991 


0.451 


1.10 


8.15 


0.892 


20 


-9.02 


9.97 


2 


16e-04 


1.12 


0.9988 


0.432 


1.12 


11.78 


0.902 


21 


-20.00 


20.00 


1.98e-04 


1.17 


0.9978 


0.395 


1.17 


26.56 


0.919 



Table 2.2: Preconditioned Richardson iterations using Theorem 2.4. 











= {ii-J + A)- 


-1 






Incomplete Cholesky 


J 




pz 


Vz 






pz 


^z 




£z 


pz 


Vz 


£z 





1 


000 


0.00 





000 


1.000 


0.00 





000 


0.643 


0.00 


0.997 


2 





510 


0.55 





751 


0.988 


0.15 





152 


0.606 


0.23 


0.997 


4 





335 


0.73 





866 


0.947 


0.33 





321 


0.554 


0.40 


0.997 


6 





226 


0.89 





926 


0.864 


0.53 





503 


0.479 


0.60 


0.998 


8 





160 


1.03 





958 


0.753 


0.72 





658 


0.391 


0.79 


0.998 


10 





122 


1.12 





973 


0.650 


0.86 





760 


0.314 


0.94 


0.998 


12 





100 


1.19 





981 


0.571 


0.96 





821 


0.256 


1.04 


0.998 


14 





087 


1.23 





985 


0.516 


1.03 





856 


0.213 


1.11 


0.998 


16 





079 


1.26 





988 


0.478 


1.07 





878 


0.182 


1.15 


0.998 


18 





073 


1.28 





989 


0.451 


1.10 





892 


0.157 


1.19 


0.997 


20 





070 


1.29 





990 


0.432 


1.12 





902 


0.138 


1.22 


0.997 
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We take ^ to be spectrally equivalent to //^Z + A, replacing the assump- 
tion (2.9) by 

m,{B;'v,v)<{{ijJ + A)v,v)<M,{B;'v,v), e V, (2.13) 
and define the associated inner product and norm, 

[v,w]:=iB;'v,w), \[v]\:=[v,vY/', (2.14) 

which now depend on z. The operator Bz{^zl + A) is then Hermitian with 
respect to [■,■], with eigenvalues in the closed interval [m^, M^]. 

In [16], the preconditioning of (2.1) by using an operator B independent 
of corresponding to /x^ = 0, was briefiy discussed, and this turned out to 
be advantageous only for small \z\. Here we shall show the following estimate 
in the present more general case for the error reduction factor with respect 
to the norm |[-]|, which is an improvement of the result in [16]. For simplicity 
we assume y = Im z > 0. 

Theorem 2.4. Consider the preconditioned equation (2.12) and the cor- 
responding iterative scheme (2.5). Let he determined as follows: With 
z = z — Hz, assume that ( = arg^ G (^tt, tt) and that B^ satisfies (2.13). Let 
(fiz — — arg az be the value in J :— {(^ — ^tt, ^tt) that maximizes the function 

J^zi^) '■— — — r-^, — — — , where A, = \z\ \\Bz\\, 

M^cos(C-(^) + A^cos(^' ^ 

and set pz — \a.z\ — i^zifz) / {'f^z cos ipz) ■ Then we have for the error reduction 
factor 

\[I - oizGzW < Sz := (1 - Vz{Vz)y'\ otz := p.e'^^^ (2.15) 
//, in addition, there is a jz ^0 such that 

Re{z[BzV,Bz{i^zI + A)v]) < -^z[BzV,vl \fv e V, (2.16) 
we define olz by choosing '{>z — — arg ocz & J to maximize the function 

U^):= , where Az:=Az-% (2.17) 

max(M2. cos(C — (/?), cos ip) \z\ 

and put Pz = l&zl = Uz{ipz)/ {mz cosifz) ■ We then have the sharper estimate 
\[I - &zGz]\ < Sz := (1 - M'fiz))'^^ - Pze-'^'. (2.18) 
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Proof. We have, for a = pe 

|[(/ - aGM? = -2R^{a[G,vM) + W?\[GM\\ (2-19) 

so, writing for brevity cq := co(</7) = cos</7 and ci := ci(</7) = cos(C — </?), and 
noting that z — e*^. 

Re (q;[Gz^;, v\)^Rea [B^(iiJ + A)v, v] + Re(az) [B^v, v] 
= p{co[B^{liJ + A)v,v] + ci\z\[B^v,v]) . 

Noting that Cq > and Ci > for (/? G J, we find, since = A^, 

\[GM\ < \[BM^zI + A)v]\ + 1^1 \[BM\ (2.21) 

< ic^^M, + cT'A,y/^{co[B,{i^J + A)v,v]+ci\z\[B,v,v]Y^\ 
so that, by (2.19) and (2.20) 

|[(/ - aG,)v]\^ < IHP - 2p{co[B,{fiJ + A)v,v] + c^\z\[B,v,v]) 

+ p\cq^M, + q'Az) {co[B,{pJ + A)v, v] + ci \z\ [B,v, v]) . 

Minimizing in p we find p — 1/{cq^Mz + Ci^h^z) — cqCi/ {c\Mz + cqA^), and 
hence 

|[(/ - aGM? < IMP - A i'^om^J + A)v,v] + c,\z\[BzV,v]). 

Here, by (2.13), 

co[Bz{pJ + A)v,v] + ci\z\[B^v,v] > co[B^{pJ + A)v,v] >m^co|[f]|^, 
and thus, remembering that p depends on ip through cq and Ci, 

l[(/ - aGM\' < IMP - 4rr^lMP = (1 - -o(^))iMP- 

Minimizing in (f over J shows the result stated. 

The first inequahty in (2.21) could be somewhat wasteful. If we assume 
that (2.16) holds, then we find, instead of (2.21), 

\[G,v]\' = \[B,{pJ + A)v]\' + \z\' \[B,v]\' + 2 Re {z[B,v, B,{pJ + A)v]) 

< \[B,{pJ + A)v]\'' + \z\^ \[B,v]\^ -2^,[B,v,v] 

< M,[B,{pJ + A)v,v] + |£|(A, - 27,/|£|)[S,^;,^;] 

< max^CQ 
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so that, by (2.19), 

I [{I-aGM? < I H I' - 2p(co[fi,(/x,/ + A)v, v] + Ci\z\ [B,v, v]) 

+ max(co ^M^, c^^K) {co[B^{iiJ + A)v, v] + Ci\z\ [B^v, v]) . 

The proof of (2.18) is now finished in the same way as that of (2.15) above. 

□ 

In the hmiting case when 2; — > and /iz — > 0, with ( — > ^tt, the method 
of analysis in Theorem 2.4 gives 



1 / TTl TTl 

«o = «o = ^ and e,^e,^^l--^l-— (2.22) 

compared to the error reduction ratio (M — m)/{M + m) ~ 1 — 2m/ M 
in (2.11). 

Applying Theorem 2.4 in the special case B^, — {nzl + A)'^, with 

ll-B^II = T-^^ , niz^Mz^l, 72 = -Re^, 

we see from Table 2.2 that, for our model problem, p^, (fz and Sz are close to 
the corresponding values in Table 2.1, but the values of Sz are worse. (The 
points Zj are the same for both tables, as are the values of Pz-) 

To better understand the condition (2.16) for general Bz, we write 

Hz Bzipzl + A) = H+ + iH; and Fz := \x - Pz\H^ - yH', (2.23) 

where the Hermitian operators are defined by 

H+:^1{Hz + H:) and H; -.^ -i\{H z - Hi) , 

and we have used * to denote the adjoint with respect to (■,■)• When Bz 
commutes with A, the operator Hz is Hermitian in V, and (2.16) follows if 
Iz < I Rez| m^;, because \i{Hz) > ruz by our assumption (2.13). This result 
is contained as the case H~ = of the following proposition. 

Proposition 2.5. Fix z — x + iy with x = x — //^ < and y > 0, and let Fz 

be the Hermitian operator defined in (2.23). Then a necessary and sufficient 
condition for (2.16) is that < 72 < Xi{Fz). 

Proof. We find that [BzV, Bz{pzl + A)v] = {v, Hzv) = {v, H+v) - i {v, H'v), 
so 

- Re{z[BzV, Bzipzl + A)v]) = \x\{v, H+v) - y {v, H;v) = {v, Fzv). 
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Table 2.3: Richardson iteration preconditioned by k V-cycles of AMG. 



j 


Pz 






7 

K 




Pz 


^z 


£z 


U 


l.UUU 


0.00 


0.643 


1 


0.000 


1 AAA 

1.000 


A A A 

0.00 


A /I O 


2 


0.517 


0.54 


0.767 


3 


0.003 


A A1 A 

0.919 


f\ A 1 

0.41 


A A £^ A 

0.464 


A 




n 7"^ 


n 87"^ 
u.o / o 


Q 

o 


u.uou 


n S1 1 


u.uo 


u.uzo 


6 


0.238 


0.88 


0.934 


2 


0.043 


0.541 


1.00 


0.869 


8 


0.166 


1.02 


0.962 


2 


0.301 


0.429 


1.13 


0.918 


10 


0.125 


1.12 


0.976 


2 


0.738 


0.341 


1.22 


0.947 


12 


0.102 


1.19 


0.983 


2 


1.351 


0.277 


1.27 


0.963 


14 


0.093 


1.23 


0.989 


1 


0.387 


0.185 


1.30 


0.983 


16 


0.083 


1.26 


0.991 


1 


1.112 


0.175 


1.32 


0.985 


18 


0.077 


1.28 


0.992 


1 


2.388 


0.169 


1.34 


0.985 


20 


0.072 


1.29 


0.992 


1 


3.426 


0.158 


1.35 


0.987 



Since is Hermitian and [BzV,v] — {v,v), it follows that 

.^^ -Re{z[B,v,B,(i^J + A)v]) _ .^^ {v,F,v) _ 
o^^vev [BzV,v] o^vev {v,v) 



□ 



In the Hermitian case, = 0, this proposition implies that (2.18) holds 
with = — 2mz\ cos(\ in (2.17). In general, since 

X,{F,)>\x\X,{H^)-y\\H;\\, 

a sufficient condition for Xi{Fz) > is that \\H~\\ < \x\y~^Xi{H^), which 
makes Hz essentially Hermitian. 

We have also the following simple consequence of Proposition 2.5. 

Corollary 2.6. // \\Hz — I\\ < 5 for some S < \x\/{\x\ +y), then (2.16) is 
satisfied with jz — \x\ — S{\x\ + y). 

Proof. Since \\H*-I\\ = we have \\H;\\ = \\l{Hz-I)-l{H*-I)\\ < 

S and \\H+ - I\\ < l\\Hz - I\\ + l\\H; - I\\ = \\Hz - I\\ < S, so it follows 
from Fz = \x\I + \x\{H+ - I) - yH' that 

{v, Fzv) > \x\ — 5\x\ llvll^ — = 7zlkll^ ^ v & V . 

Hence, Ai(F^) > 7^ > 0. □ 
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We now consider the practical application of these methods to the linear 
system (1.10). Putting Az := zM.+S so that AzW — g, the basic Richardson 
iteration (2.1) takes the form 

= (X - aM-U^)wj" + aM'^g = w?" + aM-V", 

where := g — {zM. +5)it>" denotes the nth residual. For the lumped mass 
method, we replace M. throughout by the corresponding diagonal matrix T>, 
whose inverse is trivial to compute. 

In the case of the special preconditioner — {jizl + ^)~^, we find that 
GzV — {^zl + A)~^AzV — {iJ'zM + S)~^AzV and so (2.5) takes the form 

^n+l ^ ^ a{flzM + 5)-V". 

We may write a general preconditioner in the form B^v = BzM.v, where Bz 
is Hermitian and positive-definite with respect to the standard unitary inner 
product on C^, since then {BzV,w) = {BzM.v, M.w) . In this way, 

^n+l ^^n^ aBzV''. 

The condition (2.13) is equivalent to 

mz{B-^u, u) < {{iizM + S)u, u) < Mz{B;^u, u) Vu e C^, 

which means that Xj{iiz-M.+S, B~^) belongs to the closed interval [m^, M^] for 
all j. In Tables 2.2 and 2.3, the values of fiz are the same as in Table 2.1, and 
for our computations we used best possible values = Xi{fizAi + S,B~^) 
and Mz = Xn(PzM +S,B-^). Note also that = Xn{M,B-^) = 

Xn{Bz, M.^^) because BzV = Xv is equivalent to M.v = XB^^v and to 
Bz{Mv) ^ XM-\Mv). 

To apply Proposition 2.5, we introduce Hermitian matrices 

Hj := i^zMBzM + liSBzM + MBzS) and n; := i^iSBzM - MBzS), 

so that Hfv = MHfv for all v e C^, and then put J^z '■= \x — iJiz\'H'^ — y'H~ 
so that FzV = M^^TzV- In this way, Xi{Fz) = Xi{J^z,M). 

Table 2.2 also shows the values of pz, fz and using Bz = {Cz^^)~^ for 
an incomplete Cholesky factorization Cz^^ ~ Pzl + A, computed using [10]. 
Although better than than no preconditioning, the error reduction factors 
are still too close to unity for the method to be of practical use. We can 
compare the values when 2; = to the optimal ones given by (2.11). In 
our case, Xi{BoS) = 0.0102 and Xn{BoS) = 1.55, so a = 1.28, k = 152.0 
and {k — 1)/{k + 1) — 0.987, compared to the values ao = cxq = 0.643 and 
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£o = £o = 0.997 given by (2.22). For j > 1, we found that Ai(F^) < at 
z = Zj, so we could not apply the second estimate (2.18) of Theorem 2.4. 

To find a better preconditioner, consider any symmetric, linear iterative 
process for the equation [Hz-M + S)v = g, of the form 

v^+^ = +Bz{g- ifizM + S)v^) , with 8^ = 6^. (2.24) 

Performing k steps of this iteration defines another linear iterative process, 

v^+i^ = + Bz,k{g - ifXzM + S)v^), (2.25) 

and the relation between B^ — B^^i and B^^k nia-y be seen from the error 
reduction operator: 

I - BzAp^zM + <S) = (X - BzipizM + S))\ 

It follows that Bjf, — Bz,k, so the /c-step process is also symmetric. The 
1-step process converges if and only if a{X — Bz{lJ>z-M. + <S)) C [—Qz-Qz] 
for some Qz < ^, because Bz{^z-M + S) \s symmetric with respect to the 
inner product {{fizAi + S)v,w); cf. Bramble [3, page 4]. In this case, the 
eigenvalues of BzifJ^z-M +S) lie in the interval [1 — QzA + Qz], or equivalently, 

{l-gz){iiJzM + S)-'v,v) < (BzV.v) < {l + gz){ifizM + S)-'v,v) 

for all V e C^, showing that Bz is positive-definite. In the same way, the 
eigenvalues of Bz,k he in the interval [(1 — Qz)^, (1 + Qz)''] and Bz,k is positive- 
definite. Thus, any symmetric and convergent linear iterative process yields 
a suitable preconditioner Bz^k, and moreover the hypothesis of Corollary 2.6 
will be satisfied for k sufficiently large, because Hz — Bz^kil-i'z-M. + S) ^ I 
a.s k ^ oo. 

Table 2.3 shows the results obtained when one step of the linear itera- 
tion (2.24) corresponds to a single V-cycle of a symmetric, algebraic multi- 
grid (AMG) solver [1], and thus (2.25) corresponds to k V-cycles. For each 
quadrature point Zj, the value of k shown is the smallest for which Xi{Fz) > 0, 
allowing application of Proposition 2.5. 

The need to compute ruz and Mz, and ideally also Xi{Fz), to determine a 
good choice of the acceleration parameter a, means that Richardson iteration 
is less attractive in practice than the Krylov methods of the next section, 
which do not suffer from this drawback, and also exhibit faster convergence. 

3 Conjugate gradient method 

Once again, assume that A is a positive definite Hermitian operator in a 
finite-dimensional complex inner product space V, and consider the equation 

AzW = g, where Az :— zl + A, z — x + iy, arg^; e (— tt, tt). (3.1) 
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Given wq, a preliminary guess for the solution w, we define the residual ro := 
g — A^Wq and the associated Krylov subspace of order n > 1, 

Vn := span{ro, A^ro, A^^Vq} = span{ro, Atq, A""Vo}, 

with Vq := {0}. Note that Vn depends on z through tq. The exact soloution 
of (3.1) satisfies 

{A,w,if) = {g,ip), WipeV. (3.2) 

As in the classical conjugate gradient method, we define the approximate 
solution Wn — Wq + Vn, with Vn G Vn, by Galcrkin's method, or 

{A^Wn,(p) = {g,(p), e K, (3.3) 

and find that Vn — Wn — wq satisfies 

{A^Vn, ^) = iA,{wn - Wo), V?) = (g , " {A^wq, ip) = (r^ V?), V99 e Vn. 

The solution of (3.3) is therefore unique, because if ro = we have 

{A^Vn,Vn) = zWVnW^ + {AVn,Vn) = 0, 

which implies f„ = 0. Hence there also exists a solution of the finite dimen- 
sional problem (3.3). The error Cn := w„ — w satisfies 

{A,en, ifi) = 0, e Vn. (3.4) 

To study the convergence of Wn, we introduce the norm 

\\\v\f ■= \z\\\v\\^ + {Av,v), (3.5) 

and note the following lemma. 

Lemma 3.1. If axgz — (p & (— tt, tt), then for all v, w & V we have 

KA^v, w)! < infill infill and \{AzV,v)\ > cos{^(j)) \\\v\\\'^. 

Proof. The first part follows at once from 

\{A;,v,w)\ < \z\ \{v,w)\ + \{Av,w)\ < \z\ \\v\\ \\w\\ + {Av,vy^'^ (Aw, w)^/^. 

Setting P := e~^'^^'^, the second part now results from 

Re{l3(A^v,v)) = Re(/3z)\\v\\'^ + Rel3{Av,v) 

> |2;|cos(|0)||t;|P + cos(|0) (Av,v) = cosd^)!!!^!!!^. 

□ 
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Using this lemma, we have the following quasi-optimality result. 

Proposition 3.2. Letw and Wn be the solutions of (3.1) and (3.3), respec- 
tively. Then, for argz = G (— vr, vr), 

Ill'Ji'n — < sec(50) inf — 

v€W0+Vn 

Proof. Lemma 3.1 and (3.4) show that, for any v E Wq + Vn, 

COS(|0)|||W„ - wlW^ < I {A^{Wn - w),Wn - \ ^ \ {Az{Wn -w),V-w)\ 
< \\\Wn — w\\\ \\\v — w\\\, 

which implies the result stated. □ 

We now proceed to generalize the classical convergence analysis of the 
CG method by allowing for the complex shift in A^. Let P„ denote the space 
of polynomials of degree at most n, with complex coefficients. 

Theorem 3.3. Letw andwn be the solutions of (3.1) and (3.3), respectively. 
U Qn £ o,nd Qn(0) = 1, then, for argz = G (— vr, vr), 

|||en||| < sec(^0) max |(5n(-s + A)| |||eo|||, where Wn — w. 

Proof. Let v := w + Qn{.Az)eQ. Since Qn{^) = 1 + AP„_i(A) with Pn-i G 
and ro = g - A^wq = -A^{wo - w) = -A^eo, we have Qn{A^)eo = 
eo — Pn-i{Az)ro. Hence v — Wq — Pn-i{Az)ro & Wo-\- Vn, and we conclude by 
Proposition 3.2 that 

cos(i0)|||e„||| < lilt' -will = |||Q„(AJeo|||. 

Since A^ is a normal operator. 



Similarly, 



||(5n(^z)eo|| < max \Qniz + X)\ \\eo\\. 



(AQn{Az)eo,Qn{A^)eQ) < max |(5„(z + A)|^ (Aeo, eo), 

Xecr{A) 



and we conclude that 

|||Qn(^Jeo||| < max \Qn{z + A)| |||eo|||, 

Xecr{A) 

which completes the proof. □ 



20 



We now introduce the Tchebyshev polynomial T„ e P„ defined by 

Tn{cos9) = cos{n9) for 6* e C, 

or, equivalently, since cos(i^) = cosh6', by T„(cosli^) = cosli(n^), and show 
the following consequence of Theorem 3.3. 

Theorem 3.4. With the above notation, we have, for e (— tt, tt), 

lllenlll < sec(i0) |T„(s^)|-^|||eo|||, where := + + 

Aat — Ai 

With arg A j + 2; G (— ^tt, |7r), j = 1,N, we may write 

1/ n , _7_____ . \/\n + Z - a/Ai + - 



VAat + Z + ^/Xl + z 

Furthermore, \r]z\ < 1 — cXj^^^ with c — c{z, Ai) > 0. 

Proof. The linear change of variables s — >■ r in the complex plane, 

T = - s)(Ai + ^) + (1 + s){Xn + z)), 

takes the real interval [—1, 1] onto the segment [Ai + ^, Ajv + parallel to 
the real axis. We note that t — when s — s^, so that, if we define 

n ( ^ ^n(g) X, + Xn + 2{z-t) 
Qn{r) := Tf-r^, with s , 

then Qn{T) G P„ and (5n(0) = 1. We thus have 

\TJs)\ 1 



max \Qn(X + z) \ — max \QniT)\ — max ,^ , ,, ,^ , 



and hence the first statement of the theorem follows by Theorem 3.3. 
Defining 9 by cosh^^ = |(e^ + e^^) = Sz and letting rjz = e^, we have 

Tnisz) = i;(cosh^) = cosh(ne) = i« + 7?;'^). 

Here rjz satisfies the quadratic equation rjz + r]~^ — 2s z, with roots 



VzM^z) = SzT V^l- = -iiV-Sz + 1 ± V-Sz - 1) ■ 

Setting rjz = 77z,-(s^), we find 
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Table 3.1: Error reduction by CG iteration. 



j 


Xj 


Vj 








\Vz\ 




u 


u.uu 


u.uu 


U.9o87 


u.uuuu 


u.uuu 


u.uuuu 


A A A 

U.UU 


2 


-0.05 


0.30 


0.9690 


0.0762 


O.OOz 


0.0762 


A AA 

0.00 


4 


n 1 o 
-O.io 


0.64 


u.y6yy 


U.i6oU 


n no 1 
O.Ooi 


n 1 1^ CO 
U.i6o2 


U.UU 


a 

U 


-0 Al. 


1 09 


Q708 


96Q8 


165 


9794 


00 

U.UU 


8 


-0.81 


1.51 


0.9711 


0.3749 


0.507 


0.3880 


0.00 


10 


-1.35 


2.12 


0.9703 


0.4605 


1.138 


0.4948 


0.00 


12 


-2.10 


2.93 


0.9686 


0.5221 


2.119 


0.5839 


0.00 


14 


-3.13 


4.01 


0.9659 


0.5646 


3.530 


0.6553 


0.00 


16 


-4.54 


5.45 


0.9622 


0.5939 


5.492 


0.7121 


0.00 


18 


-6.45 


7.38 


0.9577 


0.6143 


8.183 


0.7577 


0.00 


20 


-9.02 


9.97 


0.9523 


0.6287 


11.850 


0.7946 


0.00 




-20.00 


20.00 


0.9364 


0.6570 


26.894 


0.8628 


0.00 



and the stated formula for rj^ follows because —s^ + 1 = 2(AAr + z)/ (Aat — Ai) 
and —Sz — l = 2{Xi+ z) / {Xn — Xi) . Furthermore, writing y/—Sz ± 1 = a±+ib± 
we have a± > with the sign of 6+ the same as that of 6_. Thus, 

2^ {a+ - a.y + {b+ - b.y 

l'^"' (a+ + a_)2 + (6+ + 6_)2 ^ ' 
and to complete the proof we put '■— {Xn + z)/ (Ai -|- ^) = O(Ajv) and use 

V'^^ -I- i 1 _|_ / 

□ 

Since Iry^l < 1, it follows that |T„(s2)|~^ ~ Slr^^l", and so Theorem 3.4 
shows linear convergence with approximately this rate. When A = L^, so 
that Xn ~ c/i~^, the error bound is thus of order (1 — c/i)". The values 
of \rjz\ shown in Table 3.1 refer to the model problem from Section 5, for 
which Ai ~ 1 and Ajv ~ 4, UUU. Comparing the \riz\ with the corresponding 
values of Sz in Table 2.1 confirms the superiority of the CG method over the 
Richardson iteration (without preconditioning). 

We now seek to precondition the CG method applied to (3.1), and con- 
sider first the special preconditioner — {nzl + A)"^. We multiply (3.1) 
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by z :— {z — iiz)~^ and to write the equation in the form 

zw + B^w — zB^g, (3.6) 

in which thus z and B^ play the roles previously taken by z and A. In 
particular, the Krylov subspaces are now 

Vn = span{ro, B^tq, 5""Vo}, with tq = zB^g - {zl + Bz)wo, (3.7) 

and the iterates are defined by 

{{zl + Bz)Wn, (f) = {zB^g, if), Wip e Vn, Wn = Wq + Vn, Vn G (3.8) 

The earlier analysis remains valid, with Sz now replaced by 



Xi + Xn + 2z 



, with Xj := {fXz + XN+i-j) \ j = 1, N, 



and correspondingly for rjz. Theorem 3.4 then shows that the error reduction 
factor is bounded away from 1, independently of Xn- 

Theorem 3.5. For the CG method (3.8) applied to equation (3.6), and for 
the norm |||v||p = \z\ + {BzV,v), we have 

\\\en\\\ < sec(|0) \Tn{sz)\-^ llleolll, with T^{sz) = |(^ + ^7"), 

where 



_ Xn + z- \/Xi + z _ 

Vz-^ / / and \r)z\ < c{z,Xi,iiz) < 1- (3.9) 

yXN + z+yXi + z 

We want to discuss how to choose /i^ to minimize \r]z\ for a given z. In 
practice we are only interested in z = Zj with Re^^ > Re^g ^ —9/2 and 
g <C Atv, so the assumption I^j + Aat] > |^ + Ai| is not restrictive. We show 
the following. 

Lemma 3.6. Let z be fixed with I^ + AtvI > |^ + Ai|. Then \rjz\, defined in 
(3.9), is as small as possible for /iz > — Ai when 



Qz 

l^z = -Ai + _^ (Aat - Ai) > -Ai, where Qz : = 
-•- Qz 



z + Xi 



z + X 



N 



< 1. 
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Proof. It follows from (3.9) that 

~ ^ V{z + Xi)/{z + Ajv) - + Ai)/(/i^ + Ajv) 
^{z + \^)l{z + Aat) + ^J{ii^ + Ai)/(/i^ + Ajv) ' 

so with ^1 + + Ai)/(2; + Aat) and r := •>/ (Af^ + \\)l{[iz + Aat), we 

obtain 

Here, (^i > and we want to choose r > so that i){r) is as large as possible. 
A short calculation shows that -(/^'(r) = implies (Ci + t)^+C2 = 2(^i +r)T, or 
r2 = Thus, the maximum is attained when + Ai)/(//2 + Ajv) = Q'z, 

or equivalently when //^ = — Ai + (Ajv — Ai) q^j (1 — □ 

Note that //^ tends to + Ai| — Ai as Ajv tends to infinity. 

Table 3.1 includes some values of |rji^|, first for the optimal /x^ determined 
by Lemma 3.6, and then (in the final column) for /i^ = 0. Comparing the 
ItTzI with the corresponding values of in Table 2.1, we see that, once again, 
the CG method is always superior to the Richardson iteration, although in 
both cases the preconditioning becomes less effective with increasing j. 

We now consider a more general preconditioned form of (1.1), as in (2.4), 
where is an Hermitian positive definite operator, so that the equation 
may now be written 

G,w = := B,g, where = B,A, = zB, + B,A. (3.10) 

Note that B^ and B^A are Hermitian with respect to := (i?~^v,w). 

We now define the Krylov subspaces by 

K := span{fo, G^fo, . . . , G^^Vq}, where ro := - G^wq = B^ro, (3.11) 
and the CG iterates Wn by 

{A^Wn, (f) ^ {g,(p), y If e Vn, where Wn^wo + Vn with v„ e Vn, (3.12) 
or equivalently, 

[G^Wn, (p\ = [gz: <P\: ^(p E Vn, wherC W„ = + Vn with & Vn- 

The existence and uniqueness of w„ follow as before, and the inequalities 
in Lemma 3.1 remain valid, with ||| • ||| defined in (3.5). The error again 
satisfies an orthogonality property, 

and the following quasi-optimality result and its proof carry over verbatim. 
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Proposition 3.7. Letw and Wn be the solutions of (3.1) and (3.12). Then 
— will < sec(|0) inf_ IH^ — wlH, f or (f) — arg z G {—tt , tt) . 

VGWO+Vn 

The proof of the error bound of Theorem 3.3 does not remain vahd, in 
general, because of the presence of the operator in the definition of the 
Krylov spaces T^, 



4 Practical implementation of the conjugate 
gradient method 

We first derive an algorithm for computing the iterates Wn in the basic CG 
method (3.3) of Section 3. In doing so, we make repeated use of the following 
result. 

Lemma 4.1. If 1 < n < N — dim{V) then the residual rn — g — A^Wn for 
(3.3) satisfies 

Tn e K+i, and {rn, </?) = 0, y if ^ K- 

Iffoi^^, there exists N* < N such that r„ 7^ for < n < N* , and r„ = 
for n> N*. 

Proof. The first conclusion is trivial if rn = 0, so wc may assume 7^ 0. 
Since rn = g - A^{wo + Vn) = ro - A^Vn and A^K C K+i, we have r„ G 
Vn+i- The orthogonality property follows at once from (3.3). If = then 
Wn = u SO that, by (3.2) and (3.3), Wj — u also for j > n, and thus rj — 
for j > n. □ 

Lemma 4.1 shows, in particular, that the residuals Tq, ri, . . . , rn-i form 
an orthogonal basis for the Krylov space Vn ii n < N*. 

We introduce a second sequence of vectors Pn, for < n < N*, recursively: 
put po :— To and, if 7^ for < A; < n, put 

Pn+i rn+i + y^^/3„kPk, where /3nk _i^£!z!±il^. (4.1) 
^ {AzPk,Pk) 



k=Q 



Here, is well-defined since Pk ^ ensures {AzPk,Pk) 7^ 0. Also, since 
Po G Vi, we have p„ £ Vn+i (when defined). For real 2; > 0, the construc- 
tion in (4.1) amounts to applying the usual Gramm-Schmidt procedure to 
construct a new basis for Vn that is orthogonal with respect to the inner 



25 



product {AzV,w). For a general complex z, the sesquilinear form {AzV,w) 
is not an inner product. Even so, we may now show that, just as for the 
classical CG method, the sum over k in (4.1) collapses to include at most 
one non-zero term. 

Lemma 4.2. Assume tq ^ 0. Then 'Pn £ Ki+i is well defined by (4.1) for 
< n < N* , andpn ^ Vn, so that Vn+i = span{pQ, . . . ,Pn}- If n >! we have 
Pnk = for < k < n — 1. It follows that, recursively, for n + 1 < N*, 

Pn+i = Tn+i + PnPn, where /3n := /3n,„ = -ill^tll^lf^. (4.2) 
We also have {AzPmPk) — for Q <k <n — 1, and hence {AzPn, ^p) — for 

(peV. 

Proof. We prove the first statement by induction over n. To begin with, note 
that Po 7^ and po G Vi since po — ro. Let 1 < n < N* and assume Pk ^ ^ 
and Pk e Vk+i for < /c < n — 1, so that p„ is well-defined by (4.1). We 
cannot have Pn ^ Vn because then r„ = p„ — X]fc=o Pn-i,kPk G Vn and so 

= by Lemma 4.1, which would mean that n > N*. 

We now observe that, by Lemma 4.1, 

(A^r„+i, (p)^ {z- z){rn+i, (fi) + (r„+i, A^(^) = (r„+i, Az(p) for all (p e K+i, 

so /3nk = for < A; < n — 1, and thus (4.2) holds. We finally show the last 
statement by induction on n. For n — 1, the definition of (3q means that 

{AzPi,po) = {Az{ri + Popo),po) = {Azri,po) +/3o(Apo,Po) = 0. 

Now let 2 < n < A?"* and assume that {AzPn-i:Pk) — ior < k < n — 2. 
Then, since {A^rn, (p) — ior (p e 

{A^Pn,Pk) = {A^rn,Pk) + /3n-i{A^pn-i,Pk) = 0, for < /c < 71 - 2, 

and we also have {AzPn,Pn-i) = {Azrn,Pn-i) + Pn-i{AzPn-i,Pn-i) = 0. □ 

Using Wn and p„ we may compute Wn+i as follows, and hence Pn+i from 
(4.2). 

Proposition 4.3. IfO<n< N* , then 

11?^ 11^ 

Wn+l ^Wn + anPn, whcrC an := Y-;—^ r. 

{AzPn,Pn) 
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Proof. Since Wn+i — Wn & Ki+i we have Wn+i — Wn — + otnPn for some 
E Vn and some scalar «„. Since {A^p^, (p) — by Lemma 4.2, and using 
(3.3), we have 

{A^(p, (f) = {A^{Wn+l - Wn), p) - an{A^Pn: ^) = {9, P) " (9, P) = 0, 

implying that p — For we find, because (r„+i, r„) = 0, 
Here, since (APn,Pn-i) = 0, 

which shows the value of q;„ stated. □ 
Note that, by Proposition 4.3, 

Tn+l ^Tn- ^z(Wn+l " W„) = r„ - QlnA^p^, (4.3) 

SO that also the r„ may be computed recursively. Since A^Pn needs to be 
computed anyway to determine q;„ and /3„ this saves one application of A^. 
We remark that for real z > the scalar is real so —an{r-n+i, A^Pn) — 
{rn+i,rn-OinAzPn) = l^n+iP and = l^n+i || Vll^nlP, which is the formula 
used in the classical CG method. 

We readily show, using (4.3) and Proposition 4.3, that 

Pn+l = (1 + /3n)Pn " CtnAzPn " /3n-lPn-l, 

which is consistent with a result of Faber and Manteuffel [5, Section F]: if a 
matrix has a complete set of eigenvectors with all eigenvalues lying on a line 
segment in the complex plane, then there exists an inner product for which 
the CG iteration yields vectors p„ that satisfy such a three-term recurrence 
relation. 

The algorithm to compute Wn suggested by Lemma 4.2 and Proposition 
4.3 then goes as follows: Given a preliminary guess wq, compute Tq = g — 
AzWq, and set po = ''"o- The iterative step for Wn and Pn known is then to 
find first Wn+i from Proposition 4.3 and then, using (4.3) to determine r„+i, 
to find pn+i from (4.2). The iterations continue until, e.g., ||wn+i — Wn\\ or 
||rn+i|| is bounded by a tolerance, or, cf. Theorem 3.3, this holds for \riz\^. 

Consider using this algorithm when (3.1) is the linear system (1.10) aris- 
ing from the semidiscrete, standard Galerkin method applied to the heat 
equation (1.6). As before, we have V — C^, A — M.~^S and {v,w) — 
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{M.v^w). Thus, each apphcation of involves multiphcation by A4~^, 
however this cost is not incurred in the computation of q;„ and since 
{AzV, w) — z{M.v, w) + {Sv, w) . 

We now turn to the preconditioned CG method, and consider first the 
special preconditioner = {[JzI + ^)^^) and the method based on refor- 
mulating (3.1) as (3.6), with iterates defined by (3.7) and (3.8). The above 
analysis and the corresponding algorithm may be applied also in this case. In 
the iteration step, we now have r„ = zBzg—{zI+Bz)wn and in the computa- 
tion of an and /3n, the inner product {A^v, w) is replaced by (^{z I + Bz)v, w). 
In matrix form, (3.6) may be written 

zw + {/i^M + Sy^Mw = z{ii,M + Sy^g, 

and for the inner product we have 

z{v, w) + (B^v, w) = z{Mv, w) + {{iJ,zM + S)~^Mv, Mw). 

In particular, this method admits a three term recurrence relation, although 
the algorithm then requires the application of (//.^A^+^S)^^, which is normally 
more expensive than that oi M.~^. This drawback holds also in the case of 
the lumped mass variant of the spatial discretization, where A4 is replaced 
by a diagonal matrix V. 

Although, as noted at the end of Section 3, the error analyses of Theorems 
3.3 and 3.4 do not carry over to preconditioned equations of the form (3.10), 
we shall nevertheless proceed to consider the CG method for such equations, 
given by (3.11) and (3.12). We derive a recursive algorithm for computing 
the Wn, and in the same way as above first show the following, in which we 
again put r„ ^ g - A^Wn- 

Lemma 4.4. The preconditioned residual rn ■=gz — GzWn = Bzrn satisfies 

Tn e K+i and [r„, if] = (rn, </?) = for all if &Vn, forl<n< N. 

If fQ 7^ 0, there exists N* < N such that r„ 7^ for < n < N* , rn — for 
n > N*. 

We define the the sequence Pn recursively, cf. (4.1): if 7^ for < /c < 
n, set 

n 

Pn+1 ■= rn+1 + ^ /3nkPk: for n > 0, with Po := fo, (4.4) 

k=0 
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where {Pno, Pni, ■ ■ ■ , Pnn) is now the solution of the lower-triangular, {n + 
1) X {n + 1) linear system 

3 

^{A:,Pk:Pj)/3nk = -{Arn+l:Pj): for < j < U. (4.5) 

k=0 

The existence and uniqueness of the /3nk follows since the diagonal entries 
{AzPn-,Pn) are non-zero for n < N*, because otherwise Pn = and we would 
have Tn G Vn and thus r„ = 0. Unfortunately, in contrast to the situation 
earlier, (3nk 7^ is possible for A; < n — 1, which requires all the pj to be 
stored. Using the definition (4.5) of the /3nk, we may now show the following 
partial analogue of Lemma 4.2. 

Lemma 4.5. If Tq ^ 0, then ^ Pn E Vn+i and {AzPn,Pk) = for < k < 
n < N*. 

Proof. The argument used for Lemma 4.2 again establishes that 7^ p„ e 
Vn+i for < n < N*, and to prove the second claim we again use finite 
induction on n: Taking n = in (4.5) gives ^00 — ~iAzri,po)/{AzPo,Po) so 

(A^PuPo) = {Azih + /3ooPo),Po) = {A^ri:Po) + /3oo(APo,Po) = 0. 

Now let 1 < n < N* and assume that {AzPk,Pj) = for < j < A; < n. For 
< J" < n, 

j n 
{AzPn+l,Pj) ^ {AzTn+l^Pj) + ^l^nk{AzPk,Pj) + ^ /3nfe(^zPfe, Pj) = 0. 

k=Q k=j+l 

This completes the induction step and thus the proof of the lemma. □ 

As a consequence of Lemma 4.5, the conclusion of Proposition 4.3 remains 
valid: 

Proposition 4.6. IfO<n<N*, then 



Wn+1 ^Wn + OinPn, where 



{AzPn,Pn) {AzPn,Pn)' 



Proof. The beginning of the proof of Proposition 4.3 goes through verbatim, 
but since [r„+i,r„] ^ 0, 

OiniAzPn^rn) = [G';j(q;„p„) , r„] = [Gz{wn+i - Wn),rn] 
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Mro ^ g - AzWo 
Po = ^0 = B^Mro 
for n = to maxjterations do 
an = {Mrn,rn)/{AzPn,P„) 

Wn+l = Wn + OinPn 

MVn+l = MVn - OtnAzPn (or MVn+l ^Q- A^n) 

Tn+l = B^MVn 

if converged then 

break 
end if 

Solve Yl=o{^zPk,Pj)Pnk = -{AzVn+uPj) foi < J < 71 
Pn+1 — 'f'n+l + Z]fc=0 /^nfePfe 

end for 

Figure 4.1: Matrix version of CG method for AzW — g, preconditioned by B^. 
and, since r„ = p„ - Y2=l Pn-i,kPk, 

n-l 

(^zPn,^n) = {AzPn:Pn) - ^ l3n-l,k{^zPn: Pk) = (^2Pn,Pn)- 

fc=0 

□ 

Again, the residuals satisfy r„+i = r„ — a^A^p^, implying that the pre- 
conditioned residuals satisfy r^+i — Tn — anGzPn- Each iteration is now more 
expensive than in the algorithm proposed by Lemma 4.2 and Proposition 4.3, 
both in CPU time and memory requirements, and one may want to restart 
the iteration every m steps for some moderate choice of m. Figure 4.1 pro- 
vides a pseudocode outline of the method in its matrix formulation, where, 
as in the discussion following Corollary 2.6, we let Az — zM. + S and allow 
Bz to be any symmetric positive definite matrix. Notice that by working 
with M-Vn instead of r„, we can avoid computing the action of 

5 A model problem 

We now describe a concrete initial boundary- value problem (1.6), mentioned 
already in the numerical examples of Sections 2 and 3, and present some 
further illustrations of our results. 

For the domain Q we took the trapezium with vertices (1, 0), (0, 1), (—1, 1) 
and (—1,0), shown in Figure 5.1. The minimum eigenvalue of — on fl is 
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-1 



1 



Figure 5.1: The domain Q. 

close to 15, so we chose the diffusivity a = 1/15 to give a time scale of order 1 
for (1.6). We chose the data Uq and / so that the exact solution is 

u{x, y,t) = {1 + x){l - X - y) sin(7r|/) (1 + 2t)e~*, 

and used continuous, piecewise linear finite elements on a quasi-uniform, un- 
structured triangulation Th of fl, generated by the program Gmsh [8]. The 
dimension of the finite element space Vh was N — 2663, and the maximum 
element diameter was h — 0.035. The extremal eigenvalues of the opera- 
tor A = M-^S were Ai = 1.01380 and Ajv = 4006.79. 

Table 5.1 shows the (discrete) L2-norm of the error in Uq^h{t) at four val- 
ues of t, for three choices of g, as well as the norm of the solution itself. We 
see that once q is about 20, the 0(/i^) error from the spatial discretization 
dominates the 0(e~^/^°^^) error from the time discretization; cf. (1.9). (In- 
terestingly, the lumped mass approximation, in which we replace the mass 
matrix by a diagonal matrix D, gave slightly more accurate results, with 
the added bonus of more favourable extremal eigenvalues: Ai = 1.01248 and 
Aat = 1387.22.) 

Figure 5.2 shows the convergence history of the CG method (without 
preconditioning) when z = Zj, for j = 15 and q = 20. Here, e„ is the 
solver error, that is, the difference between the nth CG iterate and the exact 
solution of the discrete problem (as computed using a direct solver [4]). As 
well as the L2 error ||e„|| and the error |||e„||| in the norm (3.5), we show the 
theoretical bound of Theorem 3.4, which is pessimistic but with roughly the 
correct error reduction factor. 
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Table 5.1: Discretization error \\Uq^h{t) — u(t)\\h- 



t 


g = 10 


g = 20 


g = 30 


Mt)\\h 


0.25 


1.3436e-02 


4.3778e-04 


4.1747e-04 


0.4452 


0.50 


6.1232e-04 


1.6260e-04 


1.7541e-04 


0.4623 


1.00 


2.2024e-04 


2.1088e-04 


2.1114e-04 


0.4206 


2.00 


1.9403e-04 


1.9411e-04 


1.9411e-04 


0.2579 
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Figure 5.2: Convergence history of CG method (no preconditioning). 
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Tabic 5.2: Iteration counts at different quadrature points. 





Richardson 


CG 




j 


INV 


AMG(3) 




INV 


IC 


AMG(l) 






U 


1 


5 


250 


1 


52 


7 


1.14e+00 


3.18e-06 


2 


i 


n 

y 


ZZ / 


r 



48 


n 
i 


1 1 0^ 1 nn 
l.ioe+UU 


o.Ube-Ub 


A 


10 


15 


235 


6 


50 


8 


1.03e+00 


2.84e-06 


6 


15 


25 


242 


7 


51 


9 


7.67e-01 


2.78e-06 


8 


24 


42 


234 


8 


50 


10 


4.39e-01 


3.03e-06 


10 


39 


56 


219 


9 


46 


11 


2.21e-01 


3.86e-06 


12 


49 


57 


184 


10 


40 


11 


1.19e-01 


6.086-06 


14 


48 


45 


149 


9 


32 


10 


7.41e-02 


1.276-05 


16 


44 


37 


98 


8 


22 


9 


5.116-02 


3.836-05 


18 


32 


26 


34 


5 


11 


5 


3.69e-02 


1.916-04 


20 


8 


7 


10 


2 


3 


2 


2.71e-02 


1.876-03 



In Table 5.2 we show iteration counts at alternate quadrature points for 
several versions of the Richardson and CG iterations. In the column headings, 
INV refers to the the special preconditioner Bz = {iJizM. + S)~^ , AMG(/c) 
refers to the algebraic multigrid preconditioner [1] with k V-cycles, and IC 
refers to the incomplete Cholesky preconditioner [10]. The first CG column 
shows the results using no preconditioner. As the acceleration parameter, 
we chose a = from Theorem 2.3 in the case of the INV preconditioner, 
and a — dtz from Theorem 2.4 for AMG(3). For both sets of Richardson 
iterations, we chose //^ as in Tables 2.1-2.3, to minimize £z from Theorem 2.3. 
Likewise, for all of the preconditioned CG iterations we chose the optimal 
value of /iz for the INV preconditioner, given in Lemma 3.6. Except for j = 0, 
the AMG(l) preconditioner for CG is almost as effective as INV, requiring 
only 11 iterations in the worst case. One could also reduce the setup cost 
for AMG by using the same //^ for several nearby quadrature points, but we 
did not investigate the tradeoff between the cost saving and a possibly slower 
convergence. 

As the stopping criterion, we used 

g-Re(zj)t 

||e„|| < ej where := 5 x for 5 = 10 and t = 1. (5.1) 

(g + l)k\Zj\ 

In this way, the estimate (1.12) ensures that the additional error in Uq^h{t) 
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due to the iterative solver is less than 5. For j' = 0, we started each iteration 
with the zero vector, but for j > 1. we used the final iterate at Zj^i as 
the starting iterate at Zj. The remaining columns of the table show the 
values of and e.j. Since the former are decreasing and the latter are 

increasing, the stopping criterion becomes easier to satisfy with increasing j, 
overcoming the deterioration in the error reduction factors of the iterative 
solvers, seen in Tables 2.1, 2.2 and 3.1. 

References 

[1] W. N. Bell, L. N. Olson, and J. B. Schroder. PyAMG: Algebraic Multi- 
grid Solvers in Python v2.0, 2011. 

[2] M. Benzi and D. Bertaccini. Block preconditioning of real-valued it- 
erative algorithms for complex linear systems. IMA J. Numer. Anal., 
28:598-618, 2008. 

[3] J. H. Bramble. Multigrid Methods, volume 294 of Pitman Research Notes 
in Mathematics. Pitman, 1993. 

[4] T. A. Davis. Algorithm 832: Umfpack, an unsymmetric-pattern multi- 
frontal method. ACM Trans. Math. Software, 30:196-199, 2004. 

[5] V. Faber and T. A. Manteuffel. Orthogonal error methods. SIAM J. 
Numer. Anal, 24:170-187, 1987. 

[6] R. W. Freund. Conjugate gradient-type methods for linear systems with 
complex symmetric coefficient matrices. SIAM J. Sci. Stat. Comput., 
13:435-448, 1992. 

[7] I. P. Gavrilyuk and V. L. Makarov. Exponentially convergent algorithms 
for the operator exponential with applications to inhomogeneous prob- 
lems in banach spaces. SIAM J. Numer. Anal, 43:2144-2171, 2005. 

[8] C. Geuzaine and J.-F. Remade. Gmsh: a three-dimensional finite ele- 
ment mesh generator with built-in pre- and post-processing facilities. 

[9] Eric Jones, Travis Ohphant, Pearu Peterson, et al. SciPy: Open Source 
Scientific Tools for Python, 2001-. 

[10] M. T. Jones and P. E. Plassmann. Algorithm 740: Fortran subroutines 
to compute improved incomplete cholesky factorizations. ACM Trans. 
Math. Software, 21:18-19, 1995. 



34 



[11] W. McLean, I. H. Sloan, and V. Thomee. Time discretization via Laplace 
tranformation of an intcgro-differential equation of parabolic type. Num. 
Math., 102:497-522, 2006. 

[12] W. McLean and V. Thomee. Time discretization of an evolution equa- 
tion via Laplace transforms. IMA J. Numer. Anal, 24:439-463, 2004. 

[13] William McLean and Vidar Thomee. Maximum-norm error analysis of 
a numerical solution via Laplace transformation and quadrature of a 
fractional order evolution equation. IMA J. Numer. Anal, 30:208-230, 
2010. 

[14] WiUiam McLean and Vidar Thomee. Numerical solution via Laplace 
transforms of a fractional order evolution equation. J. Integral Equations 
AppL, 22:57-94, 2010. 

[15] D. Sheen, I. H. Sloan, and V. Thomee. A parallel method for time- 
discretization of parabolic equations based on contour integral represen- 
tation and quadrature. Math. Comp., 69:177-195, 1999. 

[16] D. Sheen, I. H. Sloan, and V. Thomee. A parallel method for time- 
discretization of parabolic equations based on Laplace transformation 
and quadrature. IMA J. Numer. Anal, 23:269-299, 2003. 

[17] V. Thomee. A high order parallel method for time discretization of 
parabolic type equations based on Laplace transformation and quadra- 
ture. Int. J. Numer. Anal Model, 2:85-96, 2005. 

[18] V. Thomee. Galerkin Finite Element Methods for Parabolic Problems. 
Springer- Verlag, Berlin, second edition, 2006. 



35 



